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ABSTRACT 

High-precision pulsar timing relies on a solar-system ephemeris in order to convert times of arrival 
(TOAs) of pulses measured at an observatory to the solar system barycenter. Any error in the 
conversion to the barycentric TOAs leads to a systematic variation in the observed timing residuals; 
specifically, an incorrect planetary mass leads to a predominantly sinusoidal variation having a period 
and phase associated with the planet's orbital motion about the Sun. By using an array of pulsars 
(PSRs J0437-4715, J1744-1134, J1857+0943, J1909-3744), the masses of the planetary systems from 
Mercury to Saturn have been determined. These masses are consistent with the best-known masses 
determined by spacecraft observations, with the mass of the Jovian system, 9.547921(2)xl(r 4 M , 
being significantly more accurate than the mass determined from the Pioneer and Voyager spacecraft, 
and consistent with but less accurate than the value from the Galileo spacecraft. While spacecraft 
are likely to produce the most accurate measurements for individual solar system bodies, the pulsar 
technique is sensitive to planetary system masses and has the potential to provide the most accurate 
values of these masses for some planets. 

Subject headings: planets and satellites: general — planets and satellites: individual (Jupiter) - 
pulsars: general 



1. INTRODUCTION 

The technique of pulsar timing can provide precise 
measurements of the rotational, astrometric, and orbital 
parameters of a pulsar by modeling the observed pulse 
times of arrival (TOAs) . The basic timing analysis pro- 
vides a fittable parametric model of delays associated 
with variations in the Euclidean distance between the 
pulsar and the Earth (resulting from Earth's orbital mo- 
tion, the proper motion of the pulsar, and its binary mo- 
tion), dispersive delays in the interstellar medium, and 
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general relativistic time dilation of clocks in the observa- 
tory and pulsar frames and al ong the propagation path 
(see, e.g. JEdwards et al~ll2006l ). The largest variable de- 
lay term is the so-called Roemer delay: the modulation 
caused by the orbital motion of the Earth relative to the 
solar system barycenter (SSB). The amplitude of this de- 
lay is up to ^500 s, while pulse TOAs for many pulsars 
are measurable with an uncertainty of much less than 
1 fj,s. This delay is co mpensated using a numerical solar 
system ephemeris (e.g.,[Standish 199cf). However, the so- 
lar system ephemerides cannot be perfect and, at some 
level, will introduce systematic effects into the timing 
process. In addition to their use in pulsar timing, these 
ephemerides are used to provide guidance information for 
space missions (in fact, this was the original motivation 
for their development), and hence there is considerable 
interest in improving their accuracy. 

The measured TOAs, ti, are related to the rotational 
phase, <pi, of the pulsar at the time of emission as follows: 

A = Wi + Y + -> C 1 ) 

where 

T i = t i+ ^ + ri) - R -^. (2) 

c 

Here, Tj is the time of pulse emission, Si and are, 
respectively, the vectors from the SSB to the geocenter 
and from the geocenter to the observatory at time tj, and 
R is a unit vector from the SSB toward the pulsar. A* 
accounts for numerous oth er delays not relevan t to the 
present discussion (see, e.g. JEdwards et aT1l2006D . Equa- 
tion (fT]) expresses the rotational behavior of the pulsar 
as a Taylor series, which for most millisecond pulsars re- 
ceives significant non-stochastic contributions from only 
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the two terms shown. Equation @ relates the times 
of emission and reception, explicitly including the varia- 
tions in light-travel time resulting from the motion of the 
observatory with respect to the SSB. If the parameters of 
the timing model are perfect, then fa is always an integer. 
The differences between the observed phase and that pre- 
dicted by the timing model are referred to as the "timing 
residuals" , usually expressed in time units through divi- 
sion by v. The best-fit timing model is generally that 
which minimizes the weighted sum of the squared resid- 
uals, where the weights are the reciprocals of the squared 
measurement uncertainties in tj. 

The most commonly used solar system ephemerides 
for pulsar timing are from NASA's Jet Propulsion Lab- 
oratory (JPL). They are constructed by numerical in- 
tegration of the equations of motion and adjustment of 
the model parameters to fit data from optical astrom- 
etry, astrolabe measurements, observations of transits 
and occultations of the planets and their rings, radar 
ranging of the planets, radio astrometry of the planets 
using very long baseline interferometry, radio ranging 
and Doppler tracking of spa cecraft, and laser ranging 
of the Moon (|Stand ish 1998). These observations con- 
strain the motion of solar system bodies with respect 
to the Earth, however they do not tightly constrain the 
planetary masses. This is reflected in the fact that the 
planetary/solar mass ratios are normally held fixed in 
the fit. 

If the vector between the observatory and SSB is not 
correctly determined, then systematic timing residuals 
will be induced. For instance, if the mass of the Jovian 
system is in error, then sinusoidal timing residuals with a 
period equal to Jupiter's orbital period will be induced. 
The identification of such residuals therefore provides a 
method to limit or detect planetary mass errors in the 
solar system ephemeris. 

In this letter we use data taken as part of the interna- 
tional effort to detect gravitational waves ([Manchester 
2008|Uanssen et al J 1200a Uenet et alll2009t Irfobbs et al. 
20101 ) using an array of pulsars to constrain the masses of 
the solar system planetary systems. These data sets are 
described in Section [2] and their analysis is discussed in 
Section [31 In Section 0] the results are presented and in 
Section [3] we discuss the potential of future observations 
and the constraints on unknown solar system bodies. 

2. DATA SETS 

The pulsars used in this analysis (listed in Table [TJ 
were chosen from the sample observed as p art of the In- 
terna tional Pulsar Timing Array project ()Hobbs et al.l 
120101) . The four pulsars were selected based upon the 
precision of their measured TOAs, the magnitude of 
timing irregularities and on the length of the data set. 
The data sets for PSRs J0437-47 15, J1744-1134, and 
J1909-3744 are those published bv lVerbiest eTail (|2009h 
except for a reweighting as described below. 

For each observation of a pulsar, typically of 1 hr du- 
ration, the data are folded at the rotation period of 
the pulsar and summed to produce a single pulse pro- 
file of relatively high signal-to-noise ratio. The TOA for 
each profile was obtained by cross-correlating the pro- 
file with a high signal-to-noise ratio template and ad- 
justing the start time of the observation for the phase 
offset between the template and observed profiles. The 



TABLE 1 
The Data Sets 



Name MJD Range Years TOAs Rms Residual (/is) 

J0437-4715 50190 - 53819 9.9 2847 0.21 

J1744-1134 49729 - 54546 13.2 342 0.64 

J1857+0943 46436 - 54507 22.1 592 1.34 

J1909-3744 52618 - 54528 5.2 893 0.17 



PSRC HIVE (jHotan et al.ll2004[ ) and tempo2 (jHobbs et al.l 
12006ft software packages were used to process the data 
and to obtain timing solutions. 

The data set for PSR J1857+0943 is a combination 
of the pre viously published T OA data from the Arecibo 
telescope (K aspi et al.l I1994D in addition to new data 
from Arecibo, Parkes, and Effelsberg. This combined 
data set is over 22 years long. Even though this data set 
is nearly 10 years longer than the other data sets in our 
sample, accurate TOAs were not obtained for just over 3 
years during the upgrade of the Arecibo telescope. The 
lack of useful data means that an arbitrary offset has 
to be included between the pre- and post-upgrade data 
sets. This arbitrary offset absorbs low-frequency power 
in the residuals which reduces the sensitivity of the fit to 
low-frequency terms. 

The uncertainties for the parameters produced by the 
standard weighted least-squares fit implemented into 
tempo2 assume that the reduced \ 2 of the fit is unity. In 
most pulsar data sets the reduced y 2 of the fit is signif- 
icantly larger than one. There are a number of possible 
reasons for this, including: radio frequency interference 
causing subtle shape changes in the profile, variations in 
the interstellar propagation path, intrinsic variations in 
the pulse profile or the pulsar rotation rate, instrumen- 
tal artifacts, errors in the clocks used to timestamp the 
data, or gravitational waves. Many of these effects have a 
steep-spectrum or "red" character and manifest approxi- 
mately as low-order polynomials in the timing residuals. 
In order to improve the estimate of the TOA uncertainty 
(and therefore the uncertainty of the parameters in the 
fit), it is common practice to introduce a multiplier that 
is applied to the TOA uncertainties at fitting. This is 
usually determined by fitting a polynomial to "whiten" 
(i.e., flatten) the residuals and then calculating the mul- 
tiplier required to bring the reduced \ 2 to unity. Because 
of the ad hoc nature of this process and because we are 
searching for long-period signals (i.e., signals with peri- 
ods similar to the length of our data sets), we use an 
improved technique to whiten the data and obtain ac- 
curate timing model parameters in the presence of red 
noise and with poorly known TOA uncertainties. This 
technique is called "Cholesky whitening" and is summa- 
rized briefly in Section 3, but will be described more fully 
in an upcoming paper. 

3. ANALYSIS 

The position of the SSB in a Euclidean frame can be 
written as a sum over all solar system bodies (including 
the Sun), where Mj is the mass of the body and bj the 
vector position of the body (where the i subscripts used 
in Equations (1) and (2) have been dropped for clarity): 

j 
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where Mt = J2j Mj ■ An erroneous set of masses Mj = 
Mj — Sj leads to an erroneous estimate of the barycenter 
vector 

^-^-5>A, (4) 

i 



where it has been assumed that 5j <C Mj and, conse- 
quently, that Y2j Mj ~ J2j Mj ■ ^ we take the origin of 
the reference frame to be at the SSB, then fog = and 



-s = -b' B 



^ 3 M T 



The error in the model time of emission is then 
(*' -s)-R 



T — T = ■ 



1 



C 



*i(6»--H), 



(5) 

(6) 
(7) 



that is, we approximate the effect of a change in the 
mass of a planet as a relocation of the SSB along the 
vector from the original SSB to that planet. The tempo2 
software package has been modified to include the right 
side of Equation as additional terms in the model. 
The modified timing model obtains bj from a specified 
version of the JPL se ries of solar system ephemerides 
(in this work, DE421. IfbTkner et al.lT2009H . The model 
parameters Sj measure the difference between the best- 
fit masses and the values assumed by the chosen solar 
system ephemeris. Indices j of 1 - 9 refer to the planets 
(and Pluto) in ascending order of mean distance from 
the Sun (note that bg = s). Examples of the induced 
timing residuals resulting from an increase in the Jovian 
system mass of 5 x 10~ 10 M Q are given in Figure [1] for 
PSRs J0437-4715 and J1857+0943 (ignoring any effects 
caused by the fitting procedures). 

In order to deal with all our data sets, we adapted 
tempo2 to fit multiple pulsars simultaneously. Fitting 
for the pulsar specific parameters is based solely on the 
TOAs of that pulsar, whereas pulsar- independent param- 
eters, such as a planetary mass, are fitted globally over 
all data sets. This procedure reduces the impact of tim- 
ing noise in individual pulsars on the derived values for 
the global parameters. 

For data sets whose post-fit residuals have a reduced x 2 
value close to 1.0, it is possible simply to fit for the plan- 
etary system mass. To determine realistic uncertainties 
for the TOAs, we selected short sections of data (~ 30 
days long depending upon the sampling) for each pulsar, 
observatory, and back-end instrument used. Weighted 
fits to each of these short sections of data gave indepen- 
dent estimates of the correction factors which were sub- 
sequently averaged and applied to the data set for that 
combination. This procedure avoided contamination of 
the correction factors by non-white noise in the data sets. 
In the presence of non-white noise, standard fitting pro- 
cedures lead to biased paramete r estimates and unde res- 
timated uncertainties (see, e.g.. iVerbiest et aLll2008j ). 

The red noise in each da ta set was modeled fo llow- 
ing the method outlined by IVerbiest et al.l (120081). For 
PSR J 0437-4715 the model developed by IVerbiest et~aT1 
(2008) was used, while for PSRs J1744-1134, 
J1857+0943, and J1909-3744, the red noise was fitted 




-5000 -4000 -3000 -2000 -1000 1000 2000 3000 4000 5000 
MJD from 50472 




-2000 



-1500 -1000 



-500 500 
MJD from 52005 



1000 1500 



2000 



Fig. 1.— Timing residuals for PSRs J1857+0943 and J0437-4715 
using the DE421 ephemeris plotted with the line indicating the 
timing signature generated by an increase in the mass of Jupiter 

of 5 x icr 10 M Q . 

by a power-law function, P oc /~ Q , with exponents of 
a = 0.9, 0.7, and 0.55, respectively. These noise models 
provide two methods to determine the parameter uncer- 
tainties. First, conservative estimates of the parameter 
uncertainties were obta ined using a Monte-C arlo simu- 
lation as described by IVerbiest et al.l (J2008). Second, 
we implemented a new technique that both whitens the 
residuals and modifies the function being fitted before 
obtaining the parameter values and uncertainties using 
a Cholesky factorization of the data covariance matrix. 
This procedure, known as "Cholesky whitening" , will be 
fully described in a forthcoming paper. 

To test our analysis technique, a new ephemeris was 
created that had identical parameters to the DE421 
ephemeris, except for a small decrease in the mass of 
Jupiter by 7 x 10~ n M . The effect of this change was 
investigated by simulating TOAs that are predicted ex- 
actly by a given timing model and the DE421 ephemeris. 
These simulated data were then analyzed using tempo2 
with the modified ephemeris. The resulting pre-fit resid- 
uals show the expected sinusoid at Jupiter's orbital pe- 
riod together with an annual term of about half the am- 
plitude of the Jupiter term. Changing the mass of Jupiter 
has many secondary effects in the modified ephemeris. 
These include a slight variation in the Astronomical Unit 
which leads to the annual sinusoid. This small effect is 
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TABLE 2 

Planetary system masses 



System 


Best-Known Mass (M ) 


Ref. 


This Work (M Q ) 


Sj/aj 


Mercury 


1.66013(7)xl0- 7 


1 


1.6584(17)xl0- 7 


1.02 


Venus 


2.44783824(4) xlCT 6 


2 


2.44783(17) xl0~ 6 


0.05 


Mars 


3.2271560(2) xlCT 7 


3 


3.226(2)xl0~ 7 


0.58 


Jupiter 


9.54791898(16) xl(T 4 


4 


9.547921(2) xl0~ 4 


1.01 


Saturn 


2.85885670(8) xlO" 4 


5 


2.858872(8) xlO" 4 


1.91 



References. 



Anderson ct al. (1987); 



(3HKonopliv et al.l (120061 ); (4) IJacobson et al.l 
(12009) . 



(2) Konop lTTeraLI (IT 999); 
(2000); (5) Jacobson ct al. 



likely to be undetectable in real data and in any case 
would be absorbed as an offset in the position of the 
pulsar by ~ 0.1 mas. The tempo2 fitting correctly re- 
covered the simulated offset in Jupiter's mass. 

4. RESULTS 

Using the DE421 ephemeris, we have obtained timing 
residuals for the four pulsars listed in Table Q] and fitted 
for a mass difference for each of the planetary systems 
from Mercury to Saturn. The resulting mass measure- 
ments are listed in Table [H where the 1-a uncertainties 
given in parentheses are in the last quoted digit. All 
results from this work are consistent with the best cur- 
rent measurements; the number of standard deviations 
between the masses derived in this work and the best- 
known masses are given in last column. 

The mass measurement for Mars was determined with- 
out the use of the PSR J0437-4715 data. A spectral 
analysis of the data set shows significant power in a broad 
feature around the period of the Martian orbit which 
could contaminate a fit for the narrow feature that would 
indicate an error in the mass of Mars. The simple red- 
noise model used to calculate the correct uncertainties 
is not detailed enough to account for this feature and so 
this data set was not used. 

Our current data sets are sensitive to mass differences 
of approximately 1O _1O M0, independent of the planet. 
Consequently, our most precise fractional mass determi- 
nation is for the Jovian system. We therefore check our 
result by comparing the Jovian system mass obtained us- 
ing different subsets of our data. In Figure^ we show the 
fitted mass difference compared with the value used for 
the DE421 ephemeris, 9.5479191563 x 10" 4 M Q , for each 
pulsar separately and the weighted mean. For compari- 
son, we also show the b est Jovian system mass from the 
Pioneer an d Voyager (|Campbell fc SvnnoTtl 119851 ) and 
the Galileo (jJacobson et al.ll2000[ ) spacecraft. The results 
obtained by fitting to individual pulsar data sets show 
a small scatter around the DE421 mass value with no 
pulsar showing more than a 2-er deviation. The weighted 
mean deviates from the best-known measurement by only 
1.1 a and has considerably smaller uncertainties than the 
mass determination derived from Pioneer and Voyager 
(ICampbell fc Svnnotti[l985l) . 

5. DISCUSSION 

While the result presented here for the Jovian system 
is more precise than the best measurement derived from 
the Pioneer and Voyager spacecraft by a factor of ^4, 
the result from the Galileo spacecraft is still better by a 
factor of ~ 20. For a pulsar timing array of 20 pulsars, 
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Fig. 2. — Values and uncertainties f or the mass of the Jovian sys - 
tem from the Pi oneer and Voyager ([ Campbell & Svnnott I1985T ), 
and the Galileo (Jacobson ct al. 2000) spacecraft , the pulsars in- 
dividually, and the array of pulsars. 
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Fig. 3. — Mass uncertainties for the Jovian and Saturnian system 
using simulated data from an array of 20 pulsars sampled every 14 
days, timed with an rms timing residual of 100 ns for different data 
spans. Also plotted are the current best mass measurements fo r 
Jupiter (Jacobson ct al. 2000) and Saturn (Jacobson et al. 2006). 



regularly sampled every two weeks, with white data sets 
and an rms timing residual of 100 ns, the uncertainty of 
the mass estimate decreases with increasing data span 
such that the mass uncertainty of the Galileo measure- 
ment for Jupiter would be reached in ~ 7 years of ob- 
servations; see the solid line in Figure [3] Note that this 
curve does not follow a simple power-law function be- 
cause of the fitting procedures that are undertaken when 
dealing with pulsar data sets. Figure [3] also shows that, 
with ~13 years of data, the uncertainty of the current 
Cassini measurement for Saturn is reached. 

These predictions rely on pulsar data sets remaining 
"white" over timescales of a decade or more at very 
high levels of timing precision. While this has yet to be 
demonstr ated, the indications from recent decade-long 
data sets (<Verbiest et al.ll2009[ ) are encouraging. 

Analysis of data from current and future spacecraft will 
produce improved measurements of planetary masses. 
For example, NASA's New Frontiers Mission to Jupiter, 
Juno, is expected to reach Jupiter in 2016 and orbit it 
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for more than a year. A major scientific objective of this 
mission is to probe Jupiter's gravitational field in detail, 
thereby providing a very precise mass for Jupiter. 

While spacecraft measurements are likely to continue 
to provide the most precise mass measurements for most 
of the planets, at least for the next decade, it should 
be noted that the pulsar measurements are indepen- 
dent with different assumptions and sources of uncer- 
tainty. Independent methods are particularly important 
for high-precision measurements where sources of sys- 
tematic error may not be well understood. Furthermore, 
while spacecraft such as Juno are very sensitive to the 
mass of the individual body being orbited, they tell us 
very little about the satellite masses of that body. Only 
five of the Jovian and nine of the Saturnian satellites 
are included in the system mass assumed for DE421 (R. 
A. Jacobson, private communication). Since the pulsar 
technique is sensitive to the mass of the entire planetary 
system, it can provide a measure of the mass undeter- 
mined by spacecraft observations. 

By combining the pulsar and satellite measurements, 
it will be possible to test the inverse-square relation of 
gravity and distance for Jupiter masses and distances 
between 0.1 and 5 AU. Although no deviations apart 
from known general relativistic effects are expected, it is 
important to place limits on such effects where possible. 

The pulsar timing technique is also sensitive to other 
solar system objects such as asteroids and currently un- 
known bodies, e.g., trans-Neptunian objects (TNOs). 
Measurements of anomalous period derivatives and bi- 
nary period derivatives for a number of millisecond pul- 
sars have already been used to place limits on the ac- 
celeration of the Solar System t oward nearby stars or 
unde t ected massive planet s (e.g., Zakam ska fe Tremaind 
120051: iVerbiest et al.ll2008h . Pulsar timing array experi- 
ments with a wide distribution of pulsars on the sky will 
be sensitive to the dipolar spatial dependence resulting 
from any error in the solar system ephemeris, includ- 
ing currently unknown TNOs. Any ephemeris error will 
be distinguishable from the effects of gravitational waves 
passing over the Earth as the latter have a quadrupolar 
spatial signature. Limits for unknown masses have also 
been placed by spacecraft using deviations from their 
predicted trajectories. Doppler tracking data from the 
two Pioneer spacecraft were searched for accelerations 
due to an unknown planet. The anomalous acceleration 
detected in these data, ap = (8.7 ± 1.3) x 10~ 10 m s -1 is 
attributed to non-gravitational sources ([Anderson et al.l 



Anderson, J. D., Colombo, G., Espsitio, P. B., Lau, E. L., k. 

Trager, G. B. 1987, Icarus, 71, 337 
Anderson, J. D., Laing, P. A., Lau, E. L., Liu, A. S., Nieto, 

M. M., & Turyshev, S. G. 2002, Phys. Rev. D, 65, 082004 
Campbell, J. K. & Synnott, S. P. 1985, AJ, 90, 364 
Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, 

MNRAS, 372, 1549 
Folkner, W. M. 2010, in IAU Symposium, Vol. 261, IAU 

Symposium, ed. S. A. Klioner, P. K. Seidelmann, & 

M. H. Soffel, 155-158 
Folkner, W. M., Williams, J. G., & Boggs, D. H. 2009, "The 

Planetary and Lunar Ephemeris DE 412", IPN Progress Report 

42-178, NASA Jet Propulsion Laboratory 
Hobbs, G. et al. 2010, CQGrav., 27, 084013 
Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, 

MNRAS, 369, 655 



I2002D and is not detected in planetary measurements 
(|Folkned[2010l) . 

An exciting possibility for the future is the creation of a 
solar system ephemeris that includes pulsar timing data 
in the overall fit. Such a fit would be able to determine 
the masses of the planetary systems while simultaneously 
fitting for orbital parameters. 

6. CONCLUSIONS 

We have used the four longest and most precise data 
sets taken for pulsar timing array projects to constrain 
the masses of the solar system planetary systems from 
Mercury to Saturn. In all cases, these measurements are 
consistent with the best-known measurements. For the 
Jovian system, our measurement improves on the Pio- 
neer and Voyager spacecraft measurement and is con- 
sistent with the mass derived from observations of the 
Galileo spacecraft as it orbited the planet between 1995 
and 2003. Pulsar timing has the potential to make the 
most accurate measurements of planetary system masses 
and to detect currently unknown solar system objects 
such as TNOs. In the future, pulsar timing data can 
be included in the global solutions used to derive solar 
system ephemerides, thereby improving their precision. 
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